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Abstract 



We review the mathematical formalism underlying the modelling of stochasticity in bio- 
logical systems. Beginning with a description of the system in terms of its basic constituents, 
we derive the mesoscopic equations governing the dynamics which generalise the more fa- 
miliar macroscopic equations. We apply this formalism to the analysis of two specific noise- 
induced phenomena observed in biologically-inspired models. In the first example, we show 
how the stochastic amplification of a Turing instability gives rise to spatial and temporal 
patterns which may be understood within the linear noise approximation. The second ex- 
ample concerns the spontaneous emergence of cell polarity, where we make analytic progress 
by exploiting a separation of time-scales. 



1. Introduction 

The evidence from the talks delivered at this meeting should leave little doubt that 
stochastic models are required to understand a wide range of biological phenomena. The 
models do need to be carefully formulated, with the nature of the interactions between the 
constituents clearly specified, but when this is carried out in an unambiguous fashion, their 
analysis can begin. In the great majority of cases this will be numerical, and many of the 
contributions at this meeting were concerned with the development and implementation of 
efficient algorithms for carrying this analysis out. 

In this paper we will take a different path. We will similarly take care in formulating the 
precise form of the model, stressing the need to begin from a microscopic individual-based 
model (IBM). However we will be primarily be interested in using mathematical analysis to 
understand the model. To facilitate this, the interactions between the constituents of the 



model will be represented by chemical-like reactions occurring at specified rates (Black and 



McKane, 2012). If the underlying stochastic process is Markov, which will usually be the 



case, a master equation can be written down which will describe the time evolution of the 



probability of the system being in a given state (van Kampen, 2007). 



One analytic procedure that can always be carried out on the master equation is to 
calculate the equations describing the time evolution of the averages of the state variables. 
These deterministic equations give a macroscopic description of the dynamics of the system. 
They are the other major methodology for the theoretical study of biological systems, and 
their use is exemplified by the book by Murray (2008 ). This older tradition involved both the 
study of simple, analytically tractable, models and dynamical systems theory. The former 
was concerned with the mathematical investigation of specific differential equations of few 
variables and the latter with general results on stability of attractors, topological notions, 
bifurcation theory, and so on (Wiggins, 2003). 

A parallel methodology for the mathematical analysis of the full stochastic model also 
exists, however it is much less widely appreciated than that for the corresponding determin- 
istic analysis. Stochastic differential equations (SDEs) can be derived when the number of 
constituents is large (but not infinite); the stochasticity originating from the discreteness of 
the underlying IBM. Techniques from the theory of stochastic processes, as well as general 
results from this theory, can be used to understand these equations analytically, just as in 
the deterministic case. 

Here we will describe this methodology with reference to two specific models in order 
to illustrate the basic ideas and techniques. We will begin in section 2 with the definition 
of the state variables of the model, and the specification of the interactions between them 
in terms of chemical reactions. This allows us to write down the master equation for the 
process. In section 3 we derive the macroscopic and mesoscopic equations governing the 
dynamics of the process from this master equation. We then go on to illustrate the use of 
this formalism, within the linear noise approximation (LNA), to the particular example of the 
Brusselator model in section 4. However, the LNA is not always able to capture the details 
of stochastic ordering phenomena, and in section 5 we show that, in some situations, we are 
able to go beyond the LNA. This is illustrated on the spontaneous emergence of cell polarity; 
this example being motivated by the talk given by Linda Petzold at the meeting. Finally 
we conclude with an overview of the techniques and applications that we have discussed in 
section 6. 



2. Definition of the models and the master equation 

The stochastic effects that will interest us here occur when the number of constituents, 
N, is large but finite. In this case a mesoscopic description is required: the state variables 
are assumed continuous — unlike in the microscopic IBM description — but the effect of 
the discreteness is not lost: it is manifested in the form of Gaussian distributed 'noise'. 
The mesoscopic form of the model is not obvious; we have to begin from a microscopic 
IBM and derive it as an approximation which holds when N is large. It cannot be found 
from a macroscopic description, since this consists of time evolution equations for average 
quantities, whereas the micro- or mesoscopic formulations describe the time evolution of 
the entire probability density function (pdf). But many stochastic processes have the same 
average, so there is in principle no unique way of determining the micro- or mesoscopic model 
from the macroscopic one. This will be especially clear later, when we derive the mesoscopic 
form of the equations and compare them to their macroscopic counterparts. 



To set up the IBM we first need to decide what are the variables which describe the state 
of the system. For the simplest case of a single type of constituent with no spatial, class, 
or other structure, it would be a single integer n = 1, 2, . . . , N representing the number of 
constituents in the system. In ecological models this could be the number of individuals 
in the population. There would then be a transition rate from state n to state n' caused 
by births, deaths, competition, predation, etc. This rate will be denoted by T(n'\n), with 
the initial state on the left (the other convention with the initial state on the right is also 
sometimes used). 

The probability of finding the system in state n at time t, P n {t), changes according to 
the master equation: 

= £ T (n\n') PAt) ~ £ T(n'\n) P n (t). (1) 

The first term on the right-hand side describes the rate at which P n (t) increases due to 
transitions into the state n from all other states n' while the second term on the right-hand 
side describes the rate at which P n (t) decreases due to transitions out of the state n to all 
other states n' . The net change then gives dP n (t)/dt. Although it is intuitively clear, it can 



also be derived from the Chapman-Kolmogorov equation for Markov processes (Gardiner 



2009). It gives the 'microscopic' description of the system and will be the starting point for 
deriving the meso- and macroscopic descriptions. 

All this generalises to several types of constituents with numbers £, m, . . . at a given 
time. Spatial structure can also be introduced with constituents located in a particular 
small volume j = 1,2... at a given time. For notational simplicity we will combine the 
index i labelling these small volumes with an index s which labels the types or classes 
of constituent, into one label J = {j, s}. Later, when carrying out the analysis of the 
differential equations we will separate the indices, and may also take the continuum limit 
in which volume labels become continuous. An additional comment is worth making at this 
stage. There is no agreed nomenclature for the small volumes: describing them as 'cells' is 
potentially confusing in a biochemical context and the term 'patches' is usually only used 
in an ecological context when the constituents are individuals. We will use the neutral term 
'domain' and talk about their 'volume' V, even when they are one- or two-dimensional, as 
will be the case with the models we discuss in this paper. 

If we now write ni for the number of constituents of a particular type in a particular 
domain, we can specify the state of the system through the vector of integers n = (ni, n 2 , . . .). 
The master equation is then simply Eq. ([I]) with n replaced by n: 

= £ T(n\n')PAt) ~ £ T(n'\n)P n (t). (2) 

Having decided what the fundamental constituents of the system are, and so how the 
states of the system are defined, the next step is to give the transition rates T(n\n'). This 
will define the model, and is best done through examples. We will naturally choose as 
examples the models that will be analysed later in the paper. The first is the Brusselator 



(Glansdorff and Prigo 


gine 


Altschuler et al. 


( 


2008 


)• 



The notation used to describe the chemical types and the rates in the Brusselator model 



follows that of Cross and Greenside (2009). In every domain % molecules of two species X 



and Y interact through the reactions of the Brusselator model (Glansdorff and Prigogine 



1971): 



A X 
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— >■ 



2X + K 4 



3X. 



X 



d 



(3) 



In order, these reactions describe: (i) the creation of a new X molecule, (ii) an X molecule 
spontaneously transforming into a Y molecule, (iii) two X molecules reacting with a Y, 
changing it to an X, and (iv) X molecules being removed from the system. The rates at 
which the reactions occur are denoted by a, b, c and d, and X and Yj are the molecules that 
are in domain i at the time that the reaction occurs. Each of these reactions are assumed to 
occur independently, and without memory of the previous states of the system. In addition 
to the reactions given in Eq. ([3]), migration reactions, which describe molecular diffusion 
from one domain to another, have to be specified. For every pair of neighbouring domains 
% and j, molecules of the two species X and Y may diffuse from one domain to the other 
according to 



X — > Xj, 



Y A Yj 



(4) 



The second example consists of a two-dimensional 'cell' enclosed within a perfectly circular 



membrane (Altschuler et al. , 2008). The centre of the cell is known as the 'cytoplasmic pool' 
and contains signalling molecules which we denote by C. In our implementation, the one- 
dimensional circular membrane is divided up into domains labelled by an index i; molecules 
lying within domains i are denoted by Mj. The following reactions take place: 



c 






Mi 


fcp|F 


c, 


Mi + C 


kfb 


2Mi 



(5) 

The first two reactions describe a molecule in the cytoplasmic pool attaching itself to a 
random domain on the membrane and a molecule detaching and returning to the cytoplasmic 
pool, respectively. The third reaction represents a molecule on the membrane attempting 
to recruit another molecule from the cytoplasmic pool to the same domain. As for the 
Brusselator, there is also diffusion, in this case along the membrane: 



Mi 



a 



M. 



3 i 



(6) 



for any pair of neighbouring membrane domains i and j. The reaction scheme is illustrated 
in Fig. 1 (after Altschuler et al). 




Figure 1: Illustration of the reactions of the cell polarisation model. The three reactions are (a) attachment 
of signalling molecules to the membrane, (b) the return of a membrane molecule to the cytoplasmic pool, 
(c) recruitment of new signalling molecules by those already on the membrane, and (d) diffusion along the 
membrane, which we implement as migration between neighbouring domains. 



Once the reaction scheme has been picked and laid out as in Eqs. ([3])-((6]), the transition 
rates T(n\n') can be chosen. This effectively specifies the model. When writing the transi- 
tion rates, we will only list the variables in the domains that are involved in the reaction. 

For the Brusselator, the number of X,i and Yi will be denoted by £j and m^, respec- 
tively. We then invoke mass action, that is, the rate of the transition for a given reaction is 
proportional to the product of the densities of the reactants, to write the transition rates as 



Ti(£i + l,mi\£i,mi) 
T 2 (£i - 1,777,3 + l\ti,™>i) 

T 3 (£i + l,mi- l\£i,mi) 



T 4 (£i - l,rrii\£i 



rrii 



£ ± 



(7) 



where the subscripts on the rates refer to the four reactions in Eq. ^ and where V is the 
volume of the domain. These transition rates are as in the usual Brusselator model (Boland 



et al. , 2009 Biancalani et al. , 2011 ), although it is worth noting that we have not imposed a 



fixed limit on the number of molecules permitted in a domain (in contrast with some previous 
work (Lugo and McKane, 2008 Biancalani et al. 2010); this is reflected in the fact we use 



the inverse volume of a domain, V~ l , rather than the inverse total number of molecules in 
a domain, as the expansion parameter). 

The reaction rates which describe molecular diffusion reactions in Eq. ^ are given by 



T 6 (mj — l,m, + l\mi,m Ji 



P 



zV 
rrij 

zV 



(8) 



Here z is the number of nearest neighbours of a domain and the index j denotes a nearest 
neighbour of domain i. 

To model spontaneous cell polarisation we denote the number of C and Mj molecules by 
I and rrii, respectively. The transition rates are then taken to be 

Ti(£-l,mi + l\£,mi) = k on — 
T 2 {£+l,m i -l\£,m i ) = k oS — , 

frn ■ 

T 3 (£-l, mi +l\£, mi ) = fcfb— 1 (9) 



for the reactions ^ and 



V 2 



T 4 (mi - l,mj + l\mi,mj) = a (10) 



for the reactions 

The transition rates for these two models can be substituted into Eq. ^ which can then, 
together with suitable initial and boundary conditions, be used to solve for P n (i). They can 



also be used as the basis for setting up a simulation using the Gillespie algorithm (Gillespie 



1976, 1977). Our approach in this paper will be to devise approximation schemes for the 
master equation and to check the results obtained from such schemes with simulations based 
on the Gillespie algorithm. 

We end this section by generalising the above formulation so that it applies to a general 
biochemical network and by extension to any network of interacting agents. To do this, 
suppose that there are L different constituents in the system. These could be labelled by 
their type, the domain they occupy, etc. We will denote them by Zj, I = 1, . . . , L and at 
a given time there will be nj of them, so that the state of the system can be specified by 
n = (ni, . . . , til), as described earlier in this section. We suppose that there are M reactions 
which interconvert species: 

L L 

J2 r i» Z i — ► X>/^/, /i = l,2,...M. (11) 

7=1 7=1 

Here the numbers r Ifl and p/ M (J = l,...,L;/i = l,..., M) describe respectively the popu- 
lation of the reactants and the products involved in the reaction. The reactions Eqs. @-(j6]) 
are simple examples of this general set of reactions. 

A quantity which is central to the structure of both the mesoscopic and macroscopic 
equations is the stoichiometry matrix, vj^ = — which describes how many molecules 
of species Zj are transformed due to the reaction fi. In the notation introduced above for 
the master equation, n' = n — v, where = (z/ lAt , . . . , v Lj ^) is the stoichiometric vector 
corresponding to reaction fi. Therefore the master equation ^ may be equivalently written 
as 

dP„(f) W 



dt 



[Tn(n\n - ^)P n -^(t) - T M (n + u„\n)P n (t)] . (12) 



Many models of interest involve only a handful of different reactions; in this situation, it is 



often convenient to rewrite the master equation as a sum over reactions as in Eq. (12), rather 
than over pairs of states n and n' as in Eq. ([2]). In the next section we will describe how the 
mesoscopic description of the models can be obtained from the master equation (12). These 



can then be used as the basis for the calculational schemes which we wish to implement. 

3. Derivation of the macro- and mesoscopic equations 

3.1. Macroscopic equation 

Before we derive the mesoscopic equations from the master equation, we will carry out 
the simpler task of deriving the macroscopic equations. This will be done for the case of the 
general biochemical network described in section |2j 



This is achieved by multiplying Eq. (12) by n, and summing over all possible values of 



n. After making the change of variable n — > n + v in the first summation, one finds that 

^^ = f>,<T,(n + .>)), (13) 

where the angle brackets define the expectation value: 

("•> = £(-)*'»(*)■ (14) 

n 

In the limit where both the particle numbers and the volume become large, we will 
take the state variables to be the concentration of the constituents ni/V, rather than their 
number rij. These will be assumed to have a finite limit as V — > oo. Specifically, the state 
of the system will be determined by the new variables 

yi = lim ni , where I — 1 , . . . , L . 



From Eq. (13) we have that 

M 



^ = E^Wm(2/), I = h--.,L, (15) 



where r = t/V and where 



fM = T J im (T^(n + i/^n)) 
lim T,({n) + v,\{n)) 



= lim T^Vy + u^Vy) 

V — »oo 



(16) 



In the above we have used the fact that in the macroscopic limit the probability distribution 
functions are Dirac delta functions and so, for instance, (n m ) = (n) m , for any integer m. 
The equation 

d f = Mv), (17) 

where 

M 

= ^1/^(1/), I = 1,...,L, (18) 



is the macroscopic equation corresponding to the microscopic master equation (12). It can 
be calculated from a knowledge of the stoichiometric matrix vj^ and the transition rates 
T M (n + u^n). We scaled time by a factor of V simply because the choice we made for the 
transition rates (|7|)-([l0|) were finite as V — > oo, but we could have easily incorporated an 
extra factor of V in these rates through a time rescaling. We also chose particularly simple 
forms for these transition rates in that they were all functions of the species concentration 
ni/V. More generally, they might separately be functions of nj and V, which become 
functions of the species concentration nj/V only when both the particle numbers and the 
volume become large, so that in the limit V — > oo they become functions of the macroscopic 
state variable y. 

3.2. Mesoscopic equation 

It is perhaps useful at this stage to recall precisely what is meant by the terms 'micro- 
scopic', 'mesoscopic' and 'macroscopic'. The microscopic description is the one based on 
the fundamental constituents whose reactions are described by relations such as those in 
Eqs. (|3])-([6]). The dynamics of the processes are described by the master equation or by the 
Gillespie algorithm. The macroscopic description has been derived from this microscopic 
description above: it only involves average quantities and their time evolution. In between 
these two levels of description is the 'meso'-level description, where the continuous variable 
of the macro-description is used, but where the stochastic effects due to the discrete nature 
of the individuals is retained. Some other authors include master equations at the meso-level 
leaving the micro-level for the world of atoms and molecules, but this does not seem such 
a useful assignment in the biological context in which we are working. The derivation of 
the mesoscopic equation follows similar lines to the calculation above, with the important 
difference that we do not take an average, or equivalently, do not take the limit V — > oo. 

We begin by substituting yj = rij/V directly into the master equation. Since, as discussed 
above, our transition rates are all functions of nj/V we simply replace T^(n + v^n) by 



f^iy) in the notation of Eq. (16). In addition we will denote the pdf P n (t) where n has been 
ilaced by Vy 

dP(y,t) 



replaced by Vy as P(y,t). With these changes we may write the master equation (12) as 

M 



dt 



E U{y -y)p(y -yit) - W)P(y,t) 



For V large, the steps is^/V are likely to be very small, suggesting that we may expand the 



functions P and / as Taylor series around y. Truncating at order 0(V 2 ), we arrive at 



dP(y,r) 
dr 



M 

E E ^ 

^=1 / 



[/ M (y)P(w,r)] 



(19) 



A4 



d 2 



Lf M (2/)P(2/,r)] 



where as before we have absorbed a factor of V into the rescaled time variable r = i/V. 



This is a Fokker-Planck equation which can be cast into the standard form (Risken |1989| 
Gardinerj [2009] ) 



o 2 



2V dyidyj 



[S/j(!/)P(l/,r)] s (20) 



where Ar(y) * s defined by Eq. (18) and where 



M 



Bu(y) = vi^vjMy), I, J = !,-■■, L. 



(21) 



In the Fokker-Planck equation (20), the continuous nature of the state variables indicates 
that the individual nature of the constituents has been lost. However, the stochasticity due 
to this discreteness has not: it now manifests itself through the function E>u(y). We can 
see this is the case through the presence of the factor 1/V. 

One might ask if this approach is consistent with the previous macroscopic derivation. 
As V — > oo, the Fokker-Planck equation reduces to the Liouville equation 

dP{y,r) d 



Or 



-Hdy-^My)P{y,r)] 



(22) 



One can check by direct substitution that the solution to this equation is P(y, r) = S(y(r) — 
y) where 5 is the Dirac delta function and where y(r) is the solution of the macroscopic 



system (17); see (Gardiner, 2009) for details. 

It is also natural to ask if it is useful to include higher order terms in V^ 1 . There are 
sound mathematical reasons for not going to higher order, for instance the pdf may become 
negative ( Risken| 1989). As we will see, for the problems that we are interested in here 
(and many others) very good agreement with simulations can be found by working with the 



Fokker-Planck equation (20). 



The Fokker-Planck equation (20) provides a mesoscopic description of the system but, 
like the master equation ^ from which it originated, it is an equation for a pdf. It is 
therefore quite distinct from the macroscopic equation (17), which is an equation for the 
state variables themselves. There does, however, exist an equation for the state variables 
which is completely equivalent to the Fokker-Planck equation (20) (Gardiner, 2009). This 
equation takes the form 

d Vl . , . 1 



dr 



My) + -i=jJ2gu(y)vAT)i 



(23) 



where the t]j{t) are Gaussian white noises with zero mean and correlator 



T 



(24) 



and where gu(y) is related to Bu(y) by 

Bu(y) = ^2giK(y)g.jK(y)- 



(25) 



K 



The mesoscopic equation (23) generalises the macroscopic ordinary differential equation 



(17) with the addition of noise terms ry(r) and so is a stochastic differential equation (SDE). 
As we will discuss below we need to specify that it is to be interpreted in the sense of Ito 
(Gardiner, 2009). Notice the direct relationship between this SDE and the macroscopic 



ODE: sending V — > oo in Eq. (23) immediately yields equation (17). 



It is important to point out that the matrices gu(y) which define the behaviour of the 
noise cannot be found from the macroscopic equations, and a knowledge of the microscopic 
stochastic dynamics is essential if one is to understand the effects of noise. It is not per- 
missible in this context to simply 'add noise terms' to the macroscopic equations to obtain 
a mesoscopic description, as some authors have done in the past. The only situation in 
which it is permissible to do this is if the noise is external to the system, that is, it does not 
originate from the internal dynamics of the system. 



We end this section with two general comments on the mesoscopic equation (23). The 



first is that while there are no strong restrictions on the form of Aj(y), there are on Bjj(y). 



From Eq. (21) we see that the matrix B is symmetric, but also that for any non-zero vector 
w. 



M 



^2 WiBijWj = ^2(w vf f^(y) > 0, 



(26) 



1:1 



since f^{y) > 0. Thus B is positive semi-definite (Mehta, 1989). It follows that B = gg 
for some non-singular matrix g, where T denotes transpose (Mehta, 1989). One way of 



constructing such a matrix is to note that since B is symmetric, it can be diagonalised 
by an orthogonal transformation defined through a matrix Ojj. Then since B is positive 
semi-definite, its eigenvalues are non-negative, and so 



B = OAO 1 



gg 



where g = OA 1 / 2 , 



(27) 



and where A and A 1//2 are the diagonal matrices with respectively the eigenvalues and square 
root of the eigenvalues of B as entries. We can take the positive roots of the eigenvalues 



without loss if generality, since the sign can always be absorbed in the rjj factor in Eq. ( 23 ) 
(its distribution is Gaussian and so invariant under sign changes). It should also be pointed 
out that we can go further and make an orthogonal transformation on the noise, (j = 



Y2i SijVJi an d leave Eq. (24), and so its distribution, unchanged. The transformation matrix 



S can then be used to define a new matrix Gu = gixSjK, so that the form of the 



mesoscopic equation (23) is unchanged. So while the procedure outlined above gives us a 
way of constructing gu(y) from Bu(y), it is not unique. 



The second comment relates to the statement made earlier, that Eq. (23) is to be inter- 



preted in the Ito sense. The singular nature of white noise means that in some cases SDEs 
are not uniquely defined by simply writing down the equation, but have to be supplemented 
with the additional information on how the singular nature of the process is to be interpreted 



(van Kampen, 2007 Gardiner, 2009). This happens when gu depends on the state of the 



system y; the noise is then said to be multiplicative. As we will see in the next section, this 
subtlety is not relevant within the LNA, since there the gjj is evaluated at a fixed point of 
the dynamics, and so ceases to depend on the state of the system. However, when going 
beyond the LNA, it is an important consideration. If one wishes to manipulate a multiplica- 



tive noise SDE like Eq. (23), then one must employ modified rules of calculus which take 



into account the contribution from the noise. We refer to Gardiner (2009) for details, and a 



complete discussion on the relationship between the Ito formulation and the alternatives. 

This completes our general discussion of the derivation and form of the mesoscopic equa- 
tion. In the next two sections we will apply it to the two models which we introduced in 
Section [2} In the first case we will consider use of the LNA is sufficient for our requirements, 
but in the second case we have to go beyond the LNA. 



4. Stochastic patterns in population models 

One of the reasons why deterministic react ion- diffusion systems are interesting is the 
fact that they may give rise to ordered structures, either in space or time. Many models 
displaying several different kinds of patterns have been extensively discussed in the literature 



(Murray, 2008 Cross and Greenside, 2009). However, the mathematical mechanisms which 



are responsible for the pattern formation are few and universal, and they can be conveniently 
analysed using the simplest models. Perhaps the most famous of these mechanisms was put 
forward by Turing in his pioneering study of morphogenesis (Turing, 1952), and it is now 



referred to as the 'Turing instability' ( Murray , 2008 ; Cross and Greenside , 2009 ) . It is invoked 



as a central paradigm in various areas of science to explain the emergence of steady spatial 
structures which typically look like 'stripes', 'hexagons' or 'spots' (Murray, 2008 Cross and 



Greenside, 2009) 



In this section, we will be interested in react ion- diffusion systems exhibiting a Turing 
instability which are composed of discrete entities as described in section [2] The intrinsic 
noise in the system will render it stochastic. As we shall see, by means of the LNA, one is 
able to make analytical progress and so clarify the role of demographic noise in the pattern 
formation process. We shall show that systems of this kind display 'stochastic patterns' 
in addition to the conventional 'Turing patterns'. It has been suggested that stochastic 
patterns are responsible for the robustness of the patterning observed in population systems 



(Butler and Goldenfeld, 2009, 2011; Biancalani et al. 2010) and they have been applied 
in several ecological models (Butler and Goldenfeld 2009 Bonachela et al. 2012), in the 



dynamics of hallucination (Butler et al. , 2012) and in a biological model with stochastic 



growth (Woolley et al. 2011). In a similar way, the emergence of stochastic travelling waves 



has been studied (Biancalani et al. , 2011 ), which has found application in a marine predator- 



prey system (Datta et al. , 2010). There is also an existing literature on stochastic patterning 



in arid ecosystems (Ridolfi et al. , 2011b) where the origin of the noise is extrinsic rather than 



intrinsic. 

The following analysis employs the Brusselator model to exemplify the general theory. 



This is a reaction scheme introduced by Lefever and Prigogine in the 1960s (Glansdorff and 



Prigogine, 1971) as a model of biochemical reactions which showed oscillatory behaviour. 



For our purposes, its interest lies in the fact that its spatial version is one of the simplest 
models which exhibits a Turing instability. 

Before we begin the analysis of this model we need to discuss some aspects of the notation 
we will use. As explained in section [2] the labels that we have been using so far combine a 
spatial index with an index for the type of constituent, for example J = {j, s}. This was 
done in order to reduce the clutter of indices. However in the analysis below we will need 
to separate them, since we will assume a regular lattice structure for the domains, which 
will allow us to use Fourier analysis to effectively diagonalise the spatial part of the system. 
The Fourier components of a given function will be labelled through the argument of that 
function, leaving only the index (e.g. s) labelling the type. Specifically we will choose the 
spatial structure to be a regular D-dimensional hypercubic lattice with periodic boundary 



conditions and domain length I. Following the conventions of (Chaikin and Lubensky, 2000) 
the discrete spatial Fourier transform is then defined as: 



D 



£ 

3=1 



-ilk-j 



J), with ./• / "il ■^<' lh ->f. 



(28) 



k=l 



where Q is the number of lattice points, j is a D-dimensional spatial vector and k is its 
Fourier conjugate. Note that % here is imaginary unit and not a spatial index, and although 
both j and k are D-dimensional vectors, we do not explicitly show this, in line with the 
notation adopted in section [2} 

We begin the analysis by deriving the well-known deterministic Brusselator equations 



from this microscopic description using Eqs. (17) and (18): 



dr 
dr 



a — (b + d) Ui + cu^Vi + ctAiii 
bui - cu^vi + /3Avi, 



(29) 



where Ui = £i/V and Vi = uii/V, where u is the density of chemical species X and v is 
the density of chemical species Y. The symbol A represents the discrete Laplacian operator 
Afj = (2/z) J2j>£dj (fj ~~ ff) wri ere j' G dj indicates that the domain j' is a nearest neigh- 
bour of the domain j and z is the co-ordination number of the lattice. The spatial Fourier 



transform of the discrete Laplacian operator reads (Lugo and McKane, 2008): 



2 
D 



D 



[cos (k s - 1] . 



(30) 



8=1 



It is possible to obtain a continuous spatial description, in the deterministic limit, by 
taking the limit of small domain length scale, I — > (Lugo and McKane, 2008). By doing so, 



one can recover the traditional partial differential equations for reaction-diffusion systems: 



du 
dv 



a — (b + d) u + cu 2 v + DiV 2 u, 



bu — cu 2 v + D 2 V 2 v, 



where D\ and D 2 are obtained by scaling the diffusivities a and (3 according to: 



(31) 



(32) 



However, we shall keep the space discrete in the following analysis because the theory is 
simpler to describe and it is most convenient for carrying out stochastic simulations. We 
shall also set / = 1, since this simply amounts to a choice of length scale, and this is the 
simplest choice. 



The macroscopic equations (29) are Eq. (17) for the particular case of the Brusselator 



model. To find the corresponding mesoscopic equations we need to find the particular form 



of Eq. (23) for the Brusselator. This we can do by calculating Bjj(y), defined in Eq. (21) 



but we will find that we do not need to utilise the non- linear equation (23) to take the 



fluctuations into account; it is sufficient to use only a linearised form. This is the LNA, and 
is implemented by writing 



yiit) = (yi(t)) 



(33) 



where (yi(t)) satisfies the macroscopic equation (17). Substituting Eq. (33) into Eq. rt23 



we expand in powers of 1/ \fV . The terms which are proportional to 1/ yV give an equation 
for f 7 : 

<U ' £ ^«W»£j + T,9iA(y))vAr), (34) 
J J 



dr 



where J is the Jacobian of the system. 

In many situations, including the one we are describing here, we are only interested in 
the fixed points of the macroscopic equation, in which case (y) = y* and the matrices J 
and g can be replaced by their values at the fixed point y*. The SDE (34) now involves only 
constant matrices: 

^=£>*^ + E^A)- (35) 

For the specific case of the Brusselator the index / includes the spatial index i and an 
index s — 1,2 which distinguishes between the variables u and v. If we take the spatial 



Fourier transform of Eq. (35), translational invariance implies that the matrices J* and g* 



are diagonalised in the spatial variables, and so this equation becomes 



d£ 7 (fc,r) 
dr 



2 

E 

6=1 



(36) 



6=1 



We are now in a position to discuss both the classical Turing patterns found in determinis- 
tic equations such as Eqs. (31 ) and the stochastic Turing patterns found in the corresponding 
mesoscopic equations. The homogeneous fixed point of Eqs. (31) is given by 



* a * M fO>7\ 

u = v = — , (37) 
a ac 



although from now on we will set c = d — 1, as is common in the literature (Glansdorff and 



Prigogine, 1971 ). In the deterministic case, the linear stability analysis about this fixed point 



is a special case of that carried out above in the stochastic picture, and corresponds to ignor- 



ing the noise term in Eq. (36). Therefore the small, deterministic, spatially inhomogeneous, 



perturbations £ 7 (fc, r) satisfy the equation 

d ^ 1 = J2^(k)Uk,r), (38) 

5=1 

where the Jacobian is found to be 

The eigenvalues of the Jacobian, A 7 (fc) (7 = 1,2), give information about the stability of 
the homogeneous state. In particular, the perturbations £ 7 (fc,r) grow like linear combina- 
tions of e^^ 1 , therefore if Re[A 7 (/c)] is positive for some k and some 7, then the perturbation 
will grow with time and the homogeneous state will be unstable for this value of k. Turing's 
insight was that the pattern eventually formed as a result of the perturbation is charac- 
terised by this value of k. The overall scenario is complicated by the nature of the boundary 



conditions, the presence of other attractors and the effect of the non-linearities (Cross and 



Greenside, 2009), but in the following we shall ignore these, and consider only the simplest 
case in order to understand the main concepts. 

In this most straightforward situation, the small perturbation which excites the unstable 
k-th Fourier mode will cause the concentrations u and v to develop a sinusoidal profile 
about their fixed point values characterised by the wave-number k. The pattern is steady or 
pulsating depending on whether or not the imaginary part Im[A 7 (/c)] is zero. In both cases, 
the amplitude of sinusoidal profile increases exponentially with a time-scale l/Re[A 7 (fc)], and 
so clearly the eigenvalue with the largest real part will dominate. By moving away from the 
homogeneous state the linear approximation will eventually lose its validity and the effect of 
the non-linear terms will become relevant. If the system admits no solutions which diverge, 
the growth will be damped by the non-linearities to some non-zero value, which defines the 
final amplitude of the spatial pattern. 

Typically, the interesting case occurs when a control parameter triggers the pattern 
formation by making the real part of one of the eigenvalues positive. For the Brusselator 
this is illustrated in Fig. [2j where the relevant eigenvalue of 3 , (i.e. the one which becomes 
positive) is shown for different values of the parameter b. Here b is the control parameter 
with the other free parameter a fixed and equal to 1.5. For b < b c ~ 2.34, the real part 
of both eigenvalues is negative and thus the homogeneous state is stable. This corresponds 




Figure 2: Real part (solid lines) and imaginary part (dashed line) of the most unstable eigenvalue of matrix 
J . Parameter values are a 
to b = 1.8 (black), b = b c = 



-- 1.5, a = 2.8 and /? 
2.34 (blue) and b = 



= 22.4, (|Cross and Greenside 
2.9 



(red). 



2009). 



Solid lines correspond 
The imaginary part is shown for b = 2.34 only, 



although it looks qualitatively the same for the range of b values displayed. 



to the situation where there are no patterns. The critical value of b, b C} occurs when the 
real part of one of the eigenvalues is tangent to the fc-axis. For values of b larger than b c , a 
window of unstable modes sets in, delimited by the intersections of Re [A] with the fc-axis. 
Each mode contributes to the pattern, although the wavelength which maximises Re [A] is 
the one with bigger amplitude as it grows faster than the other modes. 

As we have mentioned already no Turing patterns emerge from the deterministic equa- 
tions if both eigenvalues of J* have a negative real part for all k, since then the perturbations 
decay away. However, when the system is described as an IBM, intrinsic noise is present 
which acts as a continuous perturbation on the homogeneous state, exciting each /c-mode. 
Given that the homogeneous state is stable for all k, every excitation decays, although each 
with a different time-scale given by l/Re[A 7 (/c)]. Thus, the closer Re[A 7 (fc)] is to zero, the 
slower the relaxation of the corresponding fc-mode. The situation is illustrated in Fig. |2j 
where the value of k at which this occurs for b = 1.8 < b c is seen to be the maximum of 
the curve of Re[A 7 (fc)] versus k. The net effect is that only a window of modes around the 
maximum of Re[A 7 (fc)] is visible in the dynamics, the others having died away. This pat- 
terns have been called stochastic Turing patterns (Biancalani et al. 2010). Figure [3] shows 
the result of numerically simulating a two-dimensional Brusselator model using the Gillespie 
algorithm in a parameter regime where the homogeneous state is stable to perturbations for 
all fc. 

From the argument given above we would expect that in order to study stochastic Turing 
patterns it would be sufficient to linearise about the homogeneous state as before, but now 
to include noise. In addition, we would expect that the detailed nature of the noise would 
be unimportant, only its strength. This justifies the use of the LNA in the analysis of 
stochastic Turing patterns and implies that the arguments which we have used can be made 
quantitative through the use of Eq. (36). To do this we first take the temporal Fourier 
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Figure 3: (Colour online) Snapshot of two-dimensional stochastic Turing patterns for species X. The system 
consists of 40 x 40 domains with periodic boundary conditions. The parameters are the same as in Fig. [2j 
except that 6 = 2 and V = 500. Simulations started close to (u*,v*) and ran for t/V = 15. We have used 
warm colours for values of u > u* and cold colours for u < u* . White pixels indicate the fixed point value, 
u . 



transform of this equation to obtain 



(40) 



-iu8, 



From Eq. (|40|) we can find an expression for the power 



where $ 7 ^(/c, u) - 

spectrum of the fluctuations which is the quantity we use to analyse the patterns: 

P 7 (*, W ) = (\i 7 (k,oo)\ 2 ) =Y1 ^(k,u)Bl(k)(&)^(k,ou), (41) 



where B*(k) is obtained by Fourier transforming in space the matrix B(y) given by Eq. (25 ) 



evaluated at the fixed point. The details of this calculation can be found in Lugo and McKane 
(2008); here we simply state the final formulae — which holds for any two-species system 



which has one spatial dimension: 



B* u (k) 
B* 12 (k) = B* 12 



B 



12 



- 2u*aAk, 
(*0 = B 



12i 



B; 2 (k)=B* 22 -2v*(3A k 



(42) 



Here, the matrix B* indicates the correlation matrix of Eq. (25) calculated at the fixed point 



for the corresponding non-spatial system. For instance, in the case of the Brusselator this 
is obtained by considering only the reactions (|3j) without those of Eqs. Q, which yields: 
B* n = 2a(l + b), B* 12 = B* 21 = -2ab and B* 22 = 2ab. 



to 
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Figure 4: (Colour online) Analytical power spectrum of species X , obtained with the same parameter values 
as in Fig. [3j and from the analytical expression ( 41 1 . 



The expression (41) for the power is plotted in Fig. 4 We have also measured the nu- 



merical power spectrum via the Gillespie algorithm, and found good agreement with the 
analytical expression, confirming that the dynamics is captured within the approximation 
scheme we have used. The power spectrum shows a peak at k ^ and u = which in- 
dicates the presence of stochastic Turing patterns of length scale characterised by k. As 
shown in previous studies (Butler and Goldenfeld 2009 Biancalani et al. 2010), one can 



compute the region of parameters for which stochastic Turing patterns arise by looking at 
when P y (k, 0) has a maximum for some non-zero k. It has been found that those regions are 
greatly enlarged, making the pattern formation a much more robust mechanism. Specifically, 
stochastic patterns may appear even for equal diffusivities, a condition for which determin- 



Goldenfeld, 2011). 



istic patterns cannot occur (Butler and Goldenfeld 2009; Biancalani et al. 2010 Butler and 



Notice that unlike their deterministic analogue, stochastic patterns are not steady but 



they continuously decay whilst they are re-created by the effect of the noise (Scott et al. 



2011 Biancalani et al. 2011). In Fig. 5 this effect is shown by means of the dynamics of a 



one-dimensional Brusselator. The noisy nature of patterns makes them hard to detect by 
the naked eye, and the emergence of a length scale only becomes clear by means of a Fourier 
analysis. 

The theory we have presented here is rather general and also applies to other types of 
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Figure 5: (Colour online) Dynamics of a one-dimensional system of 50 domains run for 2* 10 3 time (r = t/V) 
units. Parameter values are the same as in Fig. [3j 



pattern instabilities. For instance, if the power spectrum showed a peak at k ^ and u^O 



the overall pattern would consist of stochastic travelling waves (Biancalani et al. , 2011). 



Finally, it should be mentioned that the amplitude of stochastic patterns scales as TJyV 
and it therefore vanishes in the limit V — > oo, in which the deterministic picture is recovered. 
Stochastic patterns arise because of the noise related to the discreteness of the populations 
and they are therefore less relevant for populations in which the number of individuals is 
macroscopic. 



5. Spontaneous Emergence of Cell Polarity 

Many important biological functions require cells to break their rotational symmetry and 
form a distinguished 'nose' and 'tail'. The emergence of this symmetry-breaking is known 
as polarisation, and the precise mechanisms responsible are not yet fully understood. A few 



years ago Altschuler et al. (2008) proposed a very simple model of polarisation in which 



signalling molecules self-recruit to the cell membrane, before diffusing away. They showed 
through simulations that, depending on the choice of model parameters, the membrane 
molecules may spontaneously aggregate. 

Further investigations have followed several different lines, including more detailed sim- 



ulations (Petzold et al, 2012) and mathematically rigorous studies (Gupta, 2012). In this 



section, we will show how the mesoscopic framework developed in previous sections may be 
used to analytically describe the spontaneous emergence of cell polarity in this model. This 
effect is stronger than those described by the LNA, and it will require a different theoretical 
approach. 



Before beginning the analysis, it is worth noticing that in this model the total number 
of molecules does not change; there is thus no need to distinguish between this and the 
volume, so we write N = V. Moreover, the variables i (giving the number of molecules in 
the cytoplasmic pool) and (the number of molecules in membrane domain %) are related 
by: £ + '£ i m i = V. 

As usual, we first explore the behaviour of the macroscopic equations. Substituting the 
transition rates (|5| and ^ into equations (17) and (16), we arrive at 



du 
dr 

dvi 
dr 



i 

u k on + (uk ib - fc off ) Vi + aAvi . 



Conservation of the total number of molecules implies that u + J2j v j = 1j so we may 
eliminate u from the above system to obtain 



dvj 
dr 



(l - ^2 v i) ( kon + Vi kfb ) ~ VikoS + a ^ Vi ■ 



(43) 



In the appropriate continuum limit, this equation agrees with that of Altschuler et al. (2008) 



As with the Brusselator, there is a homogeneous fixed point. Putting v % = v and setting 
dvi/dr = 0, we find the quadratic equation 



(1 - Qv)(k on + v fcft,) -vk oS = 0. 



(44) 



For simplicity, we consider the case in which k on w (that is, almost all membrane molecules 
exist as result of the feedback mechanism). In this limit, the homogeneous fixed point is 
given by v = v*/Q, where v* is the mean fraction of molecules on the membrane: 



k 



oir 



if fcft, > fc off 
otherwise. 



(45) 



To gain further insight, we pass to Fourier space, where Eq. (43) with k on = becomes 

dv k 



dt 



Vk 



fcfb(l — I 1 ^o) — ^ofr + oc( cos(lk) — l) 



(46) 



The Jacobian for this system is diagonal, and it is straightforward to read off the eigenvalues 
at the fixed point Vk = S Wt o I v* as 



&off - fcfb 



if k = 



a 



(cos(Zifc)-l) if k^O. 



We can conclude from this analysis that provided k^ > k Q g the homogeneous fixed point is 
non-zero and stable. This is a puzzle: the homogeneous state corresponds to the signalling 



molecules being spread uniformly around the membrane, if this state is stable, then how can 
polarisation occur? 

We postulate that the answer lies in the following observation. Notice that if the diffusion 
coefficient a is small then the modes with wave number k ^ are only marginally stable; 
their associated eigenvalues are close to zero. In this regime, a small random perturbation 
(resulting from intrinsic noise, for example) may be enough to push the system very far from 
its equilibrium state. Moreover, the stochastic dynamics in this regime cannot be understood 
within the framework of the LNA, for the simple reason that when the system has been 
pushed far from its steady state, linearisation around that state is no longer representative 
of the true dynamics. To make analytical progress, we will need to deal with the non-linearity 
of the model some other way. 

We begin by writing down the mesoscopic equations. For our purposes, the Fokker-Planck 



equation (20) is the most useful, with A and B given by 



A,(v) 



f 1 ~ V l) ( k ™ + V i fe fb) ~ V i k °« + aAv i ' 

3 

(l ~ Vk) (k on + Vi fcfb) + Vik Q 



Note that we have neglected terms of order a/V from the noise, as we are interested in 
behaviour when a is small and V large, so a/V is negligible. As with the macroscopic 
equations, this system is easier to analyse in Fourier space. We introduce the distribution 
P(v, t) of Fourier variables v k , which satisfies the Fokker-Planck equation 



dP(v,r) 
dr 



dv k 



A k (v)P(v,r) 



(47) 



1 d d 



B k u(v)P(v,T) 



where 



A k {v) 
B k y(v) 



v k 



fcfb(l — I 1 v ) — k oS + a(cos(lk) — l) 



v k +k' I 



/cfb(l - I t> ) + k, 



oir 



(48) 



Note that the Fourier modes may take complex values, and we use djdvy to denote differ- 
entiating with respect to the complex conjugate. For later convenience, we assume Q is odd 
and number the modes by k G { — (Q — l)/2, . . . , (Q — l)/2}, so that v_ k = v k . 

Our remaining analysis is informed by two observations. First, we note that the non- 
linearity in equation (48) arises only from terms involving vq. Second, in the interesting 
regime fcfb — k D s ^> a, we have that the eigenvalues of the macroscopic system satisfy Ao 
X k < and thus Vq is (comparatively) very stable near v*. This implies a separation of time- 
scales in the problem: we expect Vq to relax very quickly to a value near its equilibrium, 
whilst the other modes fluctuate stochastically on much slower time-scales. Combining these 



facts suggests the following strategy: we restrict our attention to only those trajectories in 
which v is held constant at Iv*. 

Conditioning on the value of Vq alters the structure of the noise correlation matrix for 



the other modes; details of the general formulation are given in Appendix B of Rogers 



remaining Fourier modes Vk with k 7^ 0, conditioned on vq taking the value Iv* 



et al. (2012b). The result is a Fokker-Planck equation for the distribution P( c \v,t) of the 

(49) 



&p(c). 
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d d 
dv k dv\, 



b£Uv)P(v,t) 



where 



A[ c \v) 



aVk ( cos(/fc) — l) 

VkVk' 



2k, 



off 



I Vk+k' — 



(50) 



Multiplying (]50|) by v k and integrating over all v we obtain a differential equation for the 

d 



mode averages: 



dr 



Thus, for all a > we have 



[v k ) = a(v k )(cos(lk) - l) . 



(51) 



(52) 



For the second-order moments the behaviour is not so trivial. Multiplying (50) this time by 
v k Vk' and integrating yields 



[vkVk') = %^(% + fc'} + {v k Vk>) a(cos(//c) + cos(lk') - 2) - ' ' "" 



dr x ' V 

The equilibrium values are thus {vkV k ') — > if k + k' 7^ mod fi, and 



1 k, 

V v* J 



1 + aV(l - cos(lk))v*/k, 



off 



(53) 



(54) 



This is our main result, although further analysis is required to interpret the implications 
for the behaviour of the model. 

In Altschuler et al. (2008) it was observed that simulations of a continuum version of this 
model exhibited the curious phenomenon of the membrane molecules grouping together, 
despite the macroscopic equations suggesting they should be spread uniformly around the 
membrane. We introduce a summary statistic to measure this effect. Suppose the membrane 
molecules are distributed according to an angular density field v(x), and let A denote the 
mean angular separation of two molecules: 

2 



A 



1 



d(x, y)(v(x)v(y)\ dx dy , 



(55) 



where 



d(x,y) 



if \x — y\ < n 



\x - y\ 



2tt — \x — y\ otherwise. 



(56) 



To compute A from our result (54) we pass to the continuum limit, taking f2 — > oo and using 
the modes Vk as the coefficients of the Fourier series of the membrane angular density field 
v(x). Taking an angular prescription for the membrane so that / = 2tt/Q, we renormalise 
by a factor of 1/2^1 and reverse the Fourier transform to obtain 



v[x) 



lim 



kx v k 



for x G [—7r, 7r) 



(57) 



fe=i 



The calculation proceeds thus: 

A = 



2tt 



\x\(y(x)v(0)) dx 
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(58) 



where if = (aVl 2 v */fc fj) 1 / 2 , which we assume to have a finite value in the limit / — > 0. The 
first equality above comes from the rotational invariance of the model meaning that we may 
fix y = 0, after which we employ Eq. (57) and Eq. (54) in turn. Now, for k ^ 



\x\e ikx dx 



-4k~ 2 




for odd k 
for even k. 



Also, the following infinite series ( |Bromwich 1926) will be useful: 

1 



k odd 



X 2 + k 2 2x 



71" / 7TX \ 

— tanh — I 



V 2 / 



Altogether, we have 
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(59) 



The limits are A(oo) = 7r/2, which corresponds to molecules spread uniformly around the 
membrane, and A(0) = 0, which corresponds to complete localisation. In figure[6]we compare 




Figure 6: (Colour online) Mean angular distance between membrane molecules, as tp is varied over several 
orders of magnitude. Red circles are the results of simulations, the blue line shows the result from Eq. (59). 
On the right are two snapshots from simulations, with (p corresponding to the points (a) and (b) in the main 
figure. Clearly (a) is polarised and (b) is not, as predicted by the theory. 



this theoretical prediction with the result of simulations for ip varying over several orders 
of magnitude. For small values of <p the membrane molecules cluster into a tight group, 
meaning that the cell has become polarised. As <p is increased (caused by an over- abundance 
of signalling molecules in relation to their rate of diffusion around the membrane) this effect 
is weakened, and the cell loses polarity. 

Finally, we discuss the physical meaning of the parameter (p. The average number 
of molecules on the membrane at equilibrium is Vv*, the typical lifetime of a membrane 
molecule is l/k Q $, and the continuum diffusion coefficient is al 2 . The quantity ip can thus 
be interpreted as the mean total distance travelled by all the membrane molecules during 
their lifetime: if this quantity is small, they must remain localised. 



6. Discussion and Conclusion 



The main aim of this paper has been to show that the analysis of stochastic models in 
biology need not only be numerical: a range of analytical techniques are available, just as 
for deterministic models. In fact the treatment of stochastic models may be simpler, since 
in many cases the noise can be considered to be a linear perturbation (the LNA) to the 
deterministic form of the dynamical equations. Linear equations such as these are easier to 
solve, especially when the fluctuations are around a stationary state. 

The deterministic, or macroscopic, description of the system is valid when the individual 
nature of the constituents is unimportant, for example when they are so numerous as to 
effectively be infinite in number. The ensemble average of the stochastic variables will also 
obey the same deterministic equation. The general form of this equation is Eq. (17), that 
is, iji = Aj(y), where the dot denotes a time derivative and the index / includes both 
a position label and another identifying the constituent type. The function Aj can be 
calculated from the elementary reactions of the process using Eq. (18). The mesoscopic, 
or stochastic, description of the system which uses the same variables as the macroscopic 
description, is in principle no different. The general form of this equation is (23), that 



is, yi = Aj(y) + V^^Ylj gu(y)r)j, where r/j is a Gaussian white noise with zero mean 
and unit strength. The only additional function which appears over and above that in the 
deterministic equation is gu, which can, like Aj, be calculated from the elementary reactions 
of the process using Eqs. (21 ) and (25). Although Eq. (23) is a straightforward generalisation 
of Eq. (17), it is much less well-known. 



There are several reasons for the perceived difficulty of using Eq. (23), probably the most 
important being the unfamiliarity of many biologists with the general theory of stochastic 
processes. We have tried to show in this paper that the stochastic theory which is required 
need not be more complicated than that of dynamical systems theory, which is applicable to 
equations such as (17). This is especially true if the LNA is a valid approximation for the 
particular system under study. If the multiplicative nature of the noise cannot be neglected, 
as in section [5j then care is required because of the singular nature of white noise. However, 
even in this case, a systematic theory has been developed that may be applicable in situations 



in which there is a separation of time-scales (Rogers et al. , 2012a Biancalani et al. , 2012 



Rogers et al, 2012b). 



We applied Eq. (23) to two sets of processes, one for which the LNA was applicable 
and one for which it was not. The former situation was discussed in section |4] where we 
revisited the problem of the emergence of spatial structures for systems of populations, 
in the paradigmatic example of the Brusselator model. Intrinsic fluctuations, which are 
intuitively thought of as a disturbing source, appear instead to be critical for the emergence 
of spatial order. More specifically, we showed how Turing patterns can arise for parameter 
values for which the macroscopic equations predict relaxation to a homogeneous state. We 
called these patterns 'stochastic patterns', as they are generated by the continuous action 
of noise present in the system. However, it can be argued that the amplitude of stochastic 
patterns might be so small that they can hardly be observed in a real population, given that 
the amplitude of the noise is small as well. Whilst this might be true for some systems, a 
recent study (Ridolfi et al. 2011a) has suggested that the response to a small perturbation 
in a pattern-forming system can be unexpectedly large, if the system displays a sufficient 
degree of 'non-normality'. The connection between non- normality and stochastic patterns 
is so far largely unexplored, and constitutes a possible further investigation in this line of 
research. 

In section[5]we discussed an example of a stochastic phenomenon which goes beyond what 
can be understood within the LNA. The stochastic patterns appearing in the Brusselator are 
noise-driven perturbations around the homogeneous state, having characteristic magnitude 
1/yV (and thus disappearing in the limit V — > oo). By contrast, the spontaneous emergence 
of cell polarity in the model of Altschuler et al requires the noise to have a more complex 
structure, which can lead the system to a state very far removed from the homogeneous 
fixed point of the deterministic equations. To characterise this process, it was necessary to 
study the full effect of the non-linear terms in the mesoscopic equations. To achieve this, we 
exploited the natural separation of time-scales occurring between the dynamics of the zeroth 
Fourier mode (which relaxes quickly to its equilibrium value) and the remaining degrees of 
freedom. This is a non-standard technique, however, it can be made relatively systematic 
and general, as will be outlined in a forthcoming paper. The LNA has played an important 
role in boosting the recognition of the importance of stochastic effects in the literature; we 



hope that methods employing the separation of time-scales may provide the next theoretical 
advance. 
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